function [Price, BSPrice, V] = EUImplicit( r, sigma, S0, K, T, steps, method, flag, NodeM )
%EUIMPLICIT Summary of this function goes here
%   Detailed explanation goes here
% method: 0=fully implicit (Default), 1 = Crank-Nicolson, 2 = C-N Rannacher
% flag = 0: put (Default), flag = 1: call
% NodeM: Node multiplier, default is 1

% Argument check
if nargin < 7
    method = 0; flag = 0; NodeM = 1; end
if nargin < 8
    flag = 0; NodeM = 1;end
if nargin < 9
    NodeM = 1; end

if method == 0
    theta = 0;
elseif method == 1
    theta = 1/2;
else
    theta = 1;
end
% Generate S grid
S = [0:0.1*K:0.4*K, 0.45*K:0.05*K:0.8*K, ...
    0.82*K:0.02*K:0.9*K, 0.91*K:0.01*K:1.1*K, ...
    1.12*K:0.02*K:1.2*K, 1.25*K:0.05*K:1.6*K, ...
    1.7*K:0.1*K:2*K, 2.2*K, 2.4*K, 2.8*K, ...
    3.6*K, 5*K, 7.5*K, 10*K];
% Delta Tao
dt = T/steps;
% Number of Asset Prices Grid = Q
Q = length(S);

% Now if we want to increase size of S grid
if NodeM ~= 1
    for i = 2:NodeM
        S_old = S;
        Q = length(S_old);
        S_new = zeros(1,2*Q-1);
        S_new(1:2:end) = S_old;
        S_new(2:2:end-1) = (S_old(1:end-1) + S_old(2:end))/2;
        S = S_new;
        Q = length(S);
    end
end

% Now get the alpha and beta vector
[Alpha, Beta] = alphaBeta(S, r, sigma);

% Now create the sparse matrix M
% TODO: Use better way to allocate the sparse matrix by pre defining
%       index sets and elements
M = spalloc(Q,Q,3*Q);
M(1,1) = r*dt;
for j = 2:Q-1
   M(j,j-1) = -Alpha(j-1)*dt;
   M(j,j) = (Alpha(j-1) + Beta(j-1) + r)*dt;
   M(j,j+1) = -Beta(j-1)*dt;
   % Now the RHS Matrix MR
end

%{
% Following are original method I used trying to get LHS matrix directyl
% Now setup Matrix M_hat
% We will use Sparce Matrix in Matlab
% Determine max non-zero points in (I+M) and (I-M)
% Because we have formula (I+M)Vn+1 = (I-M)Vn
% We denote matrix (I+M) = ML and (I-M) = MR
% Now from p105 in Notes, theta =1/2 for CN and 0 for Fully Implicit
nz = length(Alpha) + length(Beta) + Q;
ML = spalloc(Q, Q, nz);
MR = spalloc(Q, Q, nz);
ML(1,1) = 1 + (1-theta)*r*dt;
ML(Q,Q) = 1;
MR(1,1) = 1 - theta*r*dt;
MR(Q,Q) = 1;
for j = 2:Q-1
   ML(j,j-1) =  (1-theta)*-Alpha(j-1)*dt;
   ML(j,j) = 1 + (1-theta)*(Alpha(j-1) + Beta(j-1) + r)*dt;
   ML(j,j+1) = (1-theta)*Beta(j-1)*dt;
   % Now the RHS Matrix MR
   MR(j,j-1) = -theta*Alpha(j-1)*dt;
   MR(j,j) = 1 - theta*(Alpha(j-1) + Beta(j-1) + r)*dt;
   MR(j,j+1) = -theta*Beta(j-1)*dt;
end
%}

% Option Value at tao = T-T = 0, Default for put option
V = zeros(Q,1);
if flag == 0
    V(Q) = 0;
    V(1:end-1) = max(K - S(1:end-1),0);
else
    V(Q) = S(end);
    V(1:end-1) = max(S(1:end-1) - K, 0);
end
% Now iterate through time steps to solve Vn
% Create sparse identity matrix
I = sparse(1:Q,1:Q,1);
if theta ~= 1 % Either Fully Implicit or CN method
    for i=1:steps
        V = (I + (1-theta)*M)\(I - theta*M)*V;
    end
else % The C-N Rannacher Timestepping
    for i=1:2
        V = (I + M)\V;
    end
    for i=3:steps
        V = (I + 0.5*M)\(I - 0.5*M)*V;
    end
end
% Now interpolate through V to find price based on K
Price = interp1(S,V,K,'nearest');

% Calculate B-S price
[bsCall, bsPut] = blsprice(S0,K,r,T,sigma);
if flag == 0
    BSPrice = bsPut;
else
    BSPrice = bsCall;
end

[delta, gamma] = Greeks(S',V, method, flag);

end

